Nurse-patient relationship for multi-period home health care routing and scheduling problem

This article proposes a novel dynamic objective function in a multi-period home health care (HHC) problem, known as the nurse-patient relationship (NPR). The nurse-patient relationship score indicating the trust a patient has for his or her care worker increases when the same people meet regularly and decreases when they are apart. Managing human resources in HHC is a combination of routing and scheduling problems. Due to computational complexity of the HHC problem, a 28-day home health care problem is decomposed into daily subproblems, and solved sequentially with the tabu search. The solutions are then combined to give a solution to the original problem. For problems with less complex constraints, the NPR model can also be solved using exact methods such as CPLEX. For larger scale instances, however, the numerical results show that the NPR model can only be solved in reasonable times using our proposed tabu search approach. The solutions obtained from the NPR models are compared against those from existing models in the literature such as preference and continuity of care. Essentially, the analysis revealed that the proposed NPR models encouraged the search algorithm to assign the same care worker to visit the same patient. In addition, it had a tendency to assign a care worker on consecutive days to each patient, which is one of the key factors in promoting trust between patients and care workers. This leads to the efficacy of monitoring patient’s disease progression and treatment.


Introduction
This study aims to analyze the nurse-patient relationship in Home Health Care (HHC) problem. As a care worker can serve multiple patients in one day, the HHC problem is a combination of routing and scheduling problems. The solution is a daily schedule or a visiting timetable of care workers to meet with the patients at home over a planning horizon (e.g. week, month) [1]. Treating patients at home helps reduce the demand of hospitality infrastructure in a medical center and prevent the risk of infections during the hospital visits. The planning and scheduling become more complicated as the demand increases, especially with the ongoing COVID-19 pandemic in the aging society in which elderly population is growing rapidly.
World Health Organization has projected that by 2050, at least 25 percent of Europe and North America's population will be older than 65 [2]. The situation does not only exist in the developed countries, but also in developing countries like Thailand [3]. Population projection for Thailand indicated that more than 30.2 percent of the population will be aged 60 and above in 2035. With this statistics, the demand for home health care services will dramatically increase in the near future. The service providers will therefore need efficient decision tools to help them create visit schedules and pairing care workers with patients.
The goal of HHC is to create visiting schedules or plans effectively by considering multiple factors such as travel costs and patients' satisfaction. HHC problems can be divided into two types by planning horizon. The first type is a single period or a daily problem which is a plan for one day. The other type is a multi-period plan which is a plan for multiple days.
For a single-period model, the problem focuses on finding an integrated routing and scheduling plan. The plan must respect constraints related to patient visiting time window, minimum worker requirements, worker skill qualifications, and patient visit time window [4,5]. The objective function for a single-period model includes operation costs, travel distances, worker active times, preferences, etc.
For a multiple-period model, the plan puts a higher priority in a scheduling task by finding suitable care workers to meet patients regularly. At the same time, a routing task guarantees that the daily assigned jobs are practical without violating the time window constraints. Several objective functions for a multi-period model have been studied e.g. continuity of care, balancing workload, etc.
To the best of our knowledge, there was no existing research work on evaluating nursepatient relationship and trust as a level or function. Research in medical science usually collects relationship data from questionnaires or interviews and their analyses focus on trust development and the outcomes of trust that affect medical care [6]. Understanding the relationship between two people is, nonetheless, not straightforward. Defining such a relationship using a mathematical formula can be very complicated, if not impossible. In this work, we propose a concept of nurse-patient relationship. Our study proposes two relationship functions: linear and sigmoidal relationship functions. To reduce the computational complexity, a 28-day HHC problem is separated into subproblems which are then solved sequentially with the tabu search. The performances of the three models, which put priority respectively on the preference, continuity of care, and nurse-patient relationship, are evaluated and compared on 20 modified Cordeau instances.
Thus, the contributions of this paper are twofold. First, we develop for the first time a dynamic model for multi-period home healthcare routing and scheduling problems with nurse-patient relationships where the relationship is evolved and modeled through a linear and sigmoid function. Second, to avoid the model computational complexity associated with the NPR, we propose to solve the multi-period HHC problem by decomposing it into subproblems which are solved sequentially with the tabu search. To verify the effectiveness and efficiency of the proposed tabu search approach, we also attempted to solve subproblems with an exact method, the IBM ILOG CPLEX Optimizer [7]. With the maximum time limit set for each call to an optimizer, we found that the solution returned by CPLEX was not necessarily optimal. Furthermore, numerical results also revealed that, for problem instances with more complex constraints, CPLEX terminated the search without finding a solution. On all problem instances that CPLEX managed to solve, the tabu search could also find the same quality solutions. Moreover, tabu search is capable of solving all 20 problem instances in a reasonable time. This shows that the proposed tabu search approach, while substantially more efficient than CPLEX, is also very robust.
The rest of this paper is organized as follows. Section 2 presents literature reviews. Section 3 gives the definition of a multi-period HHC problem. Section 4 explains our methodology for solving the problem. Section 5 presents experimental results and analysis. Finally, Section 6 gives conclusions and future research directions.

Literature review
In this section we give the relevant literature review of the home health care problems. Arising in 1980s in order to improve the logistic of home health care, the HHC problem is an extension of the vehicle routing problem (VRP) combining routing and scheduling problems into one [1,8]. Making decisions by matching care workers with patients, the assigned care workers deliver care services to the patient's home. The HHC problem differs from the ordinary VRP due to the presence of some features of decision making [9]. This includes, for instance, continuity of care, temporal dependency, preferences, skill/qualification, frequency of visits, etc. Nevertheless, there are some features in VRP that should be considered in the HHC problem such as time-dependent traffic information [10,11].
Recent research related to the home health care problem has shifted their focus towards assignment factors. A fatigue factor is becoming increasingly important especially during the COVID-19 pandemic, and has been regarded as a crucial objective of the model. Amindoust et al. [12] considered nurse fatigue factors to generate a nurse schedule using a hybrid genetic algorithm. Zhuang and Vincent [13] investigated a nurse scheduling problem with a newly introduced labor law in Taiwan. The changes include a lower limit of weekly working hours, an increment of the minimum interval between two working shifts spanning over two days, and an introduction of a minimum period for inter-work rest. Guo and Bard [14] proposed a column generation-based method to find a midterm nurse schedule. The plan considers nurse preferences and overtime. Sarkar et al. [15] introduced an algorithm for a nurse scheduling problem, where patient recovery is treated as the main objective, and other factors such as nurse skills, and logistic conditions, as constraints. Huang et al. [16] proposed intelligent algorithms to create a personnel schedule to accommodate charge nurses and general nurses whose roster requirements are different under hierarchical management. Zĩzõvić et al. [17] suggested a scoring matrix to assign medical professionals to SAR-COV-2 hospitals.
There were several versions of mixed integer programming models for a single period HHC. Trautsamwieser and Hirsch [18] investigated a daily schedule for home health care planning where the problem was formulated using a mathematical model. The single period planning was generated by an exact solver and a variable neighbourhood search (VNS). Their studies showed that a large-scale problem instance cannot be solved by the exact solver within a reasonable time. However, the VNS algorithm was capable of finding feasible solutions for the large-scale problems, and obtaining optimal solutions for smaller problem instances. Laesanklang and Landa-Silva [19] proposed a decomposition method to solve a large-scale realworld workforce scheduling and routing problem. The method started by decomposing the whole problem into multiple subproblems based on geographical regions, and each subproblem was subsequently solved by the exact method. The solutions of subproblems were combined in order to form the solution of the original problem where a repair process may be needed to create a feasible solution. Nasir and Dang [20] proposed a variate neighbourhood search heuristic in order to create a daily home health care plan because a mixed integer programming was unable to solve a large size problem. The goal of their research was to optimize the existing routes while considering compatibility, time restrictions, contract duration, idle time and workload balance. Castaño and Velasco [21] investigated a daily home health care delivery. To find a visit plan while balancing workloads, a network flow-based model and an optimization solver were applied with computation time limit of 3600 seconds to solve a realworld problem in Columbia. Li et al. [22] studied a HHC scheduling and routing problem for outpatient services in the community care center considering time windows, skills, and working regulations as constraints. The problem was formulated as a mixed-integer, nonlinear, convex programming model incorporating travelling costs, waiting time, and patients' preference objective functions.
Tabu search is one of the most popular metaheuristics algorithms for the home health care problem. Cordeau et al. [23] applied a tabu search algorithm to solve a periodic vehicle routing problem. The method explores the solution space by finding new solutions from the current solution to the best solution in a subset of its neighborhood. This method applied two search operators: remove a customer node from one route and insert it into another route, and replace part of the route with the combination from the other route. Triki et al. [24] applied tabu search algorithm in the two-phase approach for a periodic home health care planning inspired by Cordeau et al. [23]. They applied tabu search to generate a weekly plan and a mixed integer programming based neighbourhood search method to create a daily plan. Gourc et al. [25] used tabu search to solve a single-period multi time-windows home health care scheduling problem. Liu et al. [26] applied a hybridization of tabu search heuristic and local search to solve a periodic home health care problem and it showed promising results. Umam et al. [27] used a hybridization of genetic algorithm and tabu search to find a solution to a flowshop scheduling problem. Schract et al. [28] also utilized a hybrid algorithm between tabu search and genetic algorithm for nurse scheduling problem. Chaieb et al. [29] applied tabu search algorithm to solve the problem of recovery room planning and scheduling during COVID-19 pandemic.
For the multiple period home health care problem, there were many researches in the literature investigating the problem [2,30]. The key difference between the multi-period problem and the single-period problem is that the patients may request multiple services spread over different days of week or month, and the care workers may work multiple days weekly. The procedure must take into account complex worker assignment and regulations [1].
One important aspect which distinguishes the multiple-period HHC from the single-period HHC is the continuity of care [9]. The continuity of care can be enforced at different levels. The extreme case, or the full continuity of care, happens when the patient must only meet with one care worker during the entire horizon. A relationship between each pair of nurse and patient, measured through the so-called preference score, is also an important factor that should be considered. Wirnitzer et al. [31] tackled a monthly home care rostering problem using a mathematical model where the problem was solved by an optimization solver on a high performance computer with one hour computation time limit. The result showed that none of their models could be solved to optimality within the allotted time, but the best solutions by the optimization solver were better than that from the manual ones. Martinez et al. [32] applied heuristic algorithm combining a MIP formulation and a greedy heuristic to solve HHC problems with continuity of care to create a one-week visiting schedule. Cissé et al. [9] provides a good review of these objective functions. Table 1  From these objectives, the preferences, the continuity of care and care worker switches reflect directly the quality of services and patient satisfaction [34]. Most home health care providers require continuity of care as a condition of their roster to avoid potential loss of information when a patient meets a new care worker. Receiving care from the same worker also ensures good service quality. There were several researches in the literature implementing continuity of care. Steeg and Schröder [31] proposed a model to minimize the number of care workers to visit a patient considering nurse-patient loyalty as the main indicator for continuity of care. Nickel et al. [2] defined patient-nurse loyalty as the number of additional care workers assigned to a patient. If the patient-nurse loyalty is equal to zero, then the patient is assigned to only one care worker which is the best possible solution. Carello and Lanzarone [36] proposed a model to minimize the cost associated with reassignments. They defined a binary variable of reassignment for each patient and each time slot, and the goal is to find a solution with the least number of different care workers per patient. The reassignment variable is equal to 1 if the assignment to the current time slot of the patient is changed from the previous time slot, and 0 otherwise. So if the patient is assigned to the same care worker at every time slot, the number of reassignment for the patient will be zero which is the best possible solution. Ikegami and Uno [40] investigated the home help staff scheduling problem with workload balance to equalize working hours amongst all care workers. Unlike the proposed nurse-patient relationship (NPR) model which evolves daily, none of these objective functions for HHC appearing in previous work are dynamic. As the continuity of care was defined as the number of different care workers assigned to one patient, the model attempts to minimize the number of workers assigned to a patient. Let us remark that it is possible to have two distinct solutions with the same continuity of care. Fig 1 presents two example schedules for one patient, each having two care workers, A and B, visiting a patient for six days. Evaluating the two schedule with continuity of care will result in the same score, as both plans assign two workers. However, in reality, Solution 1 should have a higher score because Worker B met the patient for five consecutive days. To resolve this issue, we propose a novel dynamic nurse-patient relationship function which continues to evolve and is updated daily with prior assignments. To the best of our knowledge, this is the first study using a dynamic function within the framework of home health care routing and scheduling problems.

Problem definition
A multi-day HHC makes assignments for care workers to serve patients at their home for multiple days in the planning horizon. Depending on the patient's condition, each patient may or may not need a visit everyday. It is also possible to have requests that are spread as specific day-of-the-week patterns. Care workers may also have days-off. All patient requests and care worker days-off during the planning horizon are known at the time a decision is made. Therefore, our goal is to create a multi-day plan for HHC that best matches a care worker with patients considering their restrictions and availabilities. The plan also takes into account travelling distances. Nevertheless, the priority is to provide service quality for patients through preference scores, continuity of care, or nurse-patient relationship.
Let us summarize now the conditions required for our model formulation.
1. Patients request dates and times of the visits.
2. All patient requests must be met.
3. Care workers may make multiple visits in one day.
4. Visits must be made during care worker availability. No overtime assignments are allowed.
5. Assignments must prioritize the care worker who is more familiar with the patient, preferably the one who made a frequent visit.
6. Assignments should guarantee that care workers can make a visit where the travel time must be taken into the decision process.
7. It is preferable to assign care workers with a higher preference score. However, a lesser preferred care workers can be assigned if no other care workers are available.
All these aspects will be investigated thoroughly in this work. We now propose mathematical programming models for HHC.

Sets and parameters
A multi-period home health care problem described in this work can be defined over a graph G = (V, E), where V = {0, 1, . . ., n} and E = {(i, j)|i, j 2 V, i 6 ¼ j} represent the vertices and the edges of the graph where i, j 2 {0, 1, . . ., n}, respectively. The depot, denoted by 0 is the initial location of all care workers, a vertex i 2 Vn{0} represents a patient's home, and n is the number of patients. In addition, C denotes the set of all care workers and D denotes the set of days in the problem planning horizon. Table 2 provides the list of notations used in this paper.

Variables
The MIP models proposed in this work have two sets of decision variables given in (1) and (2). Define a binary variable x cd ij to be an assignment for care worker c that travels from patient i to patient j on day d and a non-negative integer variable a cd j to be the arrival time of care worker c at patient j on day d: a cd j 2 Z 0þ ; the time care worker c arrives at patient j on day d:

Constraints
Constraint (3) guarantees the coherence of each route, i.e. if care worker c visits patient j, then the care worker must leave that visiting location. For constraint (4), a route must begin and end at the depot only once each day.
Constraint (5) guarantees that all patient requests are served where d d j is a binary parameter equal one if patient j requests a service on day d, and zero otherwise. Constraint (6) ensures that assignments on day d are assigned to care workers who are on duty on day d. Note that the on-duty parameter z cd is also a binary parameter, equal one if the care worker c is available on day d, and zero otherwise. Constraint (7) prevents time window violations. Constraint (8) assigns arrival time to visit patient j. Constraint (9) guarantees that the arrival time a cd j is within the care worker availability.

Objective function
In this work, we consider and compare three models, namely, Basic, continuity of care (CC), and nurse-patient relationship (NPR).
3.4.1 Basic model. In the Basic model, we minimize a linear combination of the total travel distances and the preference costs. The objective function is given by (12). The total travel distances are a summation of distances computed from all care workers while the total preference costs are a summation of the preference score contributed from all care workers. The higher the preference score, the lower the objective function value.

Continuity of care (CC) model.
Built upon f Basic , the continuity of care model considers the maximization of continuity of care in addition [1,31]. The equivalent minimization problem is given by (13). The term relating to the continuity of care cost simply counts the total number of care workers that serve each patient. In this case, if a patient is served by the same care worker, the objective value is lower.
The MIP model for this case requires an additional parameter which is a continuity of care weight w 3 . Also, a binary variable y c j is defined to be 1 if the care worker c has been assigned to visit patient j at least once, as shown in (14)- (16).
if care worker c has been assigned to visit patient j at least once

Nurse-patient relationship (NPR) model. (1) NPR with linear relationship function and its limitations.
To capture the acquaintance level, we introduce a variable t c;d j that increases as the care worker c regularly meets the patient j, and decreases otherwise. Based on the artificial pheromone in ant colony optimization algorithm [42], t c;d j is computed using the formula in (17), where Q is a growth constant and ρ is a decay rate of the relationship when a care worker has not met a patient for several days: where the initial relationship score is zero, i.e. t c;0 j ¼ 0; 8j 2 V; c 2 C. Let us remark that despite its simplicity, there are some theoretical limitations with the linear NPR function. First, since the linear function has no bound, the linear NPR term grows without bound as we move forward in the planning horizon, making it difficult to control the weights of objective functions. This also prevents one from carrying the NPR relationship into the next planning horizon. Also, another limitation of the linear relationship is the strict assumption of linearity of the relationship over time which may not necessarily hold true for the relationship.
(2) NPR with sigmoidal relationship function. The NPR with sigmoidal relationship framework, introduced in this section, is extensible to the next planning horizon and possesses some desired properties. Essentially, not only do sigmoid curves increase over time (with bound), they also represent the law of diminishing marginal utility which states that the marginal utility of a service declines as more of it is received by an individual.
The sigmoidal relationship score is derived as the summation of the nurse-patient relationship score of care workers to be deployed on each day. The solution whose assignments are made with a higher nurse-patient relationship score will have a lower objective function value. This is represented by the additional term of the objective function in (18).
where w 4 is a relationship weight, t c;d j quantifies the level of acquaintance, and f ðt c;d j Þ is a relationship score between care worker c and patient j on day d.
Our assumption for the nurse-patient relationship score is that the relationship level will increase when a care worker visits a patient repeatedly. In contrast, if the care worker is absent or other care workers meet the patient, the relationship with the absent care worker is deteriorating. Moreover, the relationship score should increase slowly during the first few visits, and increase more rapidly once the care worker and patient get better acquainted. Finally, the relationship value should not exceed some predefined upper bound. We thus assume that our nurse-patient relationship follows the sigmoid function or S-curved [43].
To this end, we define the relationship score f ðt c;d j Þ, with the zero lower bound (the lowest relationship score) and the upper bound of 1 (the highest relationship score) as shown in (19) where k is a slope, and b is a shifting constant: Note that for the linear relationship function, the function f ðt c;d j Þ in (18) is simply the identity map.
Let us now dive into discuss time complexities of the NPR models with sigmoid and linear functions. There is no doubt that the sigmoid function will require a longer computation time compared to the linear function. It is known however that using the floating point exponentiation, the complexity is always Oð1Þ. For this reason, there is not a significant difference between the time complexity for the two relationship functions, be it linear or sigmoid. In particular, the overall complexity of a typical step of the algorithm is dominated by the procedure used to find the solutions in the search space (Algorithm 3), therefore the difference in the computation time required for the two relationship functions can in fact be negligible. Empirical evidence that shows the computation time to solve NPR models with sigmoid and linear relationship functions are not significantly different can be found in Table 6.

Comment: Computational complexity of the multi-period HHC problem
A single-period home health care problem is classified as an NP-hard problem [2]. Thus, the multi-period home health care problem is also an NP-hard problem with additional dimensions for the planning periods. Our mathematical programming model for the multi-period problem has two sets of decision variables: a set of binary variable x cd ij which allocates visits for care worker c to visit patient j after patient i on day d, and a set of variable a cd j which captures visiting time for care worker c when visiting patient j on day d. Considering only these two sets of variables, the total number of decision variables in the objective function is jCj � jDj � jVj 2 þ jCj � jDj � jVj ¼ jCj � jDj � jVj � ðjVj þ 1Þ: Solving the problem as a whole, the search space for the binary variables is 2 jCj�jDj�jVj 2 combinations. Here, the number of combinations is referred to as the maximum number of subproblems generated during the branch-and-bound operations in conventional MIP solvers. Solving the problem considering all time periods at once also requires impractically large amounts of memory. In particular, for the multi-period HHC model, the memory required for the LP relaxation could take up to 1 GB. Obviously, the amount of memory is increased during branch-and-bound operations as the number of potential subproblems due to branching can be as many as 2 jCj�jDj�jVj 2 .
In contrast, considering the problems day-by-day concerns only 2 jCj�jVj 2 combinations per subproblem. The memory required here in this case is significantly lower than that required to solve the problem as a whole. In particular, the memory required to solve each subproblem is approximately 100 MB for the LP relaxation, where the maximum number of branch-andbound subproblems is 2 jCj�jVj 2 subproblems. It was evident that using optimization solvers may require extremely high computational resources. For example, Wirnitzer et al. [31] used high performance computing with 121 GB of Ram to solve a problem consisting of 37 nurses, 143 patients, a total of 1114 visits for 7 days with 3 shifts per day. This vast amount of resources is not available in a personal computer. Thus, we propose in the next section to tackle this by decomposing the multi-period problem into multiple daily subproblems, reducing significantly the computational requirements. While mixed integer nonlinear programming solvers exist in the literature, such algorithms do not guarantee optimality when tackling a non-convex function, and a metaheuristic such as tabu search has more popularity for solving this problem.

Methodology
In this section, we present a method to solve multi-period HHC problems. Although solving the complete set of multi-period problem as a whole guarantees finding the optimal solution, the method is very computationally expensive as discussed. To reduce the computational complexity, the method given in Algorithm 4 separates the problem into |D| subproblems where | D| is the number of days in the multi-day problem. The algorithm solves the subproblems sequentially and generates |D| sub-solutions. The solutions of the |D| subproblems are then combined to give a solution to the original problem. For each subproblem, the method generates an initial solution using a greedy algorithm, and the tabu search algorithm is then used to improve the solutions. The assignments made by the tabu search will be lead by the three models presented in Section 3.4. For the NPR model, after the sub-solution for a 1-day assignment is improved, the relationship scores are updated, following the procedure described in the previous section, before the process to solve the subsequent subproblem begins. In this way, the relationship score can be viewed as a parameter of the subproblem. The processes are repeated until all subproblems are solved. Update solution initSol d 7:

Initial solution
Update L d : remove l j from L d 8: end if 9: end for 10: end while 11: return initSol d � For the Basic model, assign the available care worker with highest preference score. � For the CC model, assign the available care worker who has previously visited patient l j and with the highest preference score.
� For the NPR model, assign the care worker with highest nurse-patient relationship score. Our method starts by generating an initial solution using the greedy algorithm. The overall procedure is presented in Algorithm 2 with some specific details for each model. For the Basic model, the algorithm assigns an available care worker with highest preference score. For the CC model, the algorithm assigns a care worker who has been assigned previously and with highest preference score. Finally, the NPR model assigns a care worker with highest nursepatient relationship score. In all three models, this assignment made in the initial solution must not violate any constraints, otherwise the next available care worker with highest preference score (for Basic and CC) or highest relationship score (for NPR) is assigned, given that the assignment does not violate any constraints. The procedure is repeated until all patients are assigned.

Tabu search
After generating an initial solution, the solution is improved by tabu search. This procedure is presented in Algorithm 3. Behaving like a local search algorithm, tabu search accepts also nonimproving solutions to escape from a local optimum trap [44]. A key feature of the tabu search is the use of temporary memory, called tabu list, which records solutions that have previously been visited. A long length of the tabu list creates diversification search which explores wider solution regions and forbids the algorithm to revisit solutions already examined. A shorter length of the tabu list intensifies the search in the area that is deemed to have good solutions. During the search, the algorithm requires search operators to generate new solutions from existing solutions [45]. The operator applied in this work is shift (1,0) move-relocating a customer from one route to another. Select two routes r 1 , r 2 in S randomly 6: Apply Shift(1,0) to r 1 , r 2 , finding the best neighborhood solution S 0 which is neither in tabu list nor violates constraints 7: S S 0 and update tabu list 8: Set i+ = 1 9: if F(S 0 )<BestF then 10: FinalSol d S 0 11: BestF = F(S 0 ) 12: i nonimprove = 0 13: else 14: Set i nonimprove + = 1 15: end if 16: end while 17: return FinalSol d

Computational experiments
In this section, we present experiments using the proposed method to solve modified benchmark problems. The algorithms were written in Python and tested on a machine with Intel Core i5 CPU @ 2.3GHz, 8GB RAM. Parameters configurations and problem instances are presented in Sections 5.1 and 5.2, respectively. The results of the three models are given in Section 5.3. The last section is devoted to an analysis of the effects of applying a relationship score.

Parameter setting
For the proposed NPR model, parameters for the benchmark problems are specified as follows. Parameters ρ and Q in the NPR function are decay rate and growth constant of the relationship. A high decay rate ρ reduces the relationship score faster if a care worker is absent. Also, a high growth constant Q increases the relationship score more when a care worker makes a frequent visit to the same patient. The parameters k and b are chosen to create a sigmoid function with range in (0, 1). While these parameters can be adjusted to accommodate different levels of relationship, Table 3 summarizes the values of the model parameters used in our numerical experiments. The values of weights assigned to each objective function are also given in the same table.
The highest priority objective function for the HHC problem is the term related to the continuity of care (w 3 of f CC ) or the nurse-patient relationship (w 4 of f NRP ). The second most important objective is the preference objective (w 2 ). Finally, the lowest priority objective is the term related to total travel distance (w 1 ). In order to prioritize the importance of each objective function, we follow the idea given in [8,19]. Here, we set the preference weight w 2 equal to the distance between the service center and the furthest location. We set w 3 in the CC model equal to 2 × w 2 so that the weighted sum objective function will prioritize the continuity of care over the preference score. Finally, the weight w 4 in the NPR model is also twice the weight of the preference score w 2 . The parameter α 0 = max{α 0,j } is the maximum distance between the depot 0 and patient locations j of each instance. Table 3. Model parameters used in numerical experiments. The parameter α 0 = max{α 0,j } is the maximum distance between the depot 0 and patient locations j of each instance.

Instances
Instances for this research are based on Cordeau instances which were originated for the classical VRPTW. Cordeau instances are stored on the website of VRP-REP [46] dedicated to the study of VRPs. The original Cordeau datasets have 20 different instances. We use the instances consisting of 25 patients and 25 care workers. The following procedure is carried out to transform the Cordeau instances into the ones appropriate for HHC.
1. Generate preference scores between every pair of care workers and patients. The range of preference score is 0 to 1. Higher scores indicate a higher preference level. A non-preferred care worker with score -1 may also be assigned to a patient if no preferred care workers are available on some days.
2. Create a multi-day instance with periodic demands. The original instance is modified to a multi-day instance with mixed periodic demands. Periodic demands mean the patients reserve cares in a pattern such as every week, every two weeks, etc. Thus, each request is distributed in the pattern of the same day-of-week.
The information of periodic instances consists of the number of patients who want to be served and the number of available care workers on each day as presented in Table 4. The first column represents the instance name, the other columns represent the day index. The two rows represent the number of the patients who want to be served on each day (#Patient), and the number of available care workers (#Worker), respectively.

Overall results
This section presents results obtained from the Basic, CC, and the two NPR relationship models. Since the results of the NPR models with linear (NPR(L)) and sigmoid (NPR(S)) relationships were not significantly different, we will discuss them together as the NPR models hereafter, and defer comments on their similarities to Subsection 5.3.1.
To compare the quality of the solutions found by our solution method, an alternative approach which employs the simplex algorithm was implemented to solve the linear NPR subproblems. In particular, the effectiveness of our approach is tested against the state-of-the-art IBM ILOG CPLEX Optimizer [7]. We shall refer to this variant as NPR(L)-Cpx. However, in order not to interrupt the flow of the presentation, we will postpone the discussion of the comparisons with NPR(L)-Cpx to Subsection 5.3.2. Table 5 presents the value of the four objectives, namely, the total travel distance (Distance), the total preference score (Preference Score), the number of different care workers (Diff Worker) and the total relationship score (Relationship Score). The comparison of the four models can be done by comparing the objective function value evaluated at the optimal solution returned by the model. Therefore, even if the model itself does not involve a particular objective function, one could still compute the objective function value at the optimal solution, see e.g. Wirnitzer et al. [31]. Here, we use the sigmoidal function as the measurement of the relationship score.
Considering the first column, we can see that the Basic model required the fewest number of trips while the NPR models required the largest number of trips. Looking at the column "Distance", overall the Basic model performs better than the other two, which came as no surprise as the basic model gives priority to minimizing the distance. As for "Preference Score", the three models achieve similar results, but the average value for the NPR is highest. As for "Diff Worker", it is also apparent that the total number of different care workers required for the CC and NPR models are less than that of the Basic model. Finally, the results in the column "Relationship Score" clearly show that the NPR models have highest relationship score for all instances while the Basic model has the lowest relationship score for all instances. Table 6 displays computational times (in seconds) to find a solution for each instance. Due to complexity of the relationship score, the NPR models required significantly more computation time than the other two models. #Worker 12 12 13 13 15 13 14 13 13 14 12 17 14 14 13 12 13 12 17 14 14 11 12 14 12 17 14   Additionally, we examine daily cumulative relationship scores of each instance. Fig 2 presents the cumulative total relationship scores of the four models over the planning horizon. We can see from the graphs that NPR models have highest cumulative relationship scores overall, followed by the CC, and Basic models, respectively.

Comparison between the NPR models with linear and sigmoid relationship functions.
Results for the two NPR models in Table 5 indicated that the NPR model with sigmoid and linear functions each achieved the best NPR score on 6 problems, and were tied on 8 problems. Thus, based on our results, the NPR models are indifferent to relationship function choice. This conclusion is later supported by a statistical test in Section 5.3.3. Taking a look at the computation times for the two NPR models, as presented in Table 6, it appears there is no significant difference. These results confirm our discussions regarding the time complexities of the two relationship functions in Section 3.4.3. It is worth reiterate that despite having nearly the same performance on this benchmark suite, practically the sigmoidal NPR is more flexible than the linear one due to its extensibility to the next time horizon.

Comparing the tabu search-based with the CPLEX-based approaches.
The objective of this section is to determine the quality of the relationship scores found with the tabu search heuristic. When solving the subproblems, the alternative method is based on the simplex algorithm, CPLEX. As we found no statistically significant difference between the linear and sigmoid NPR models, we will only consider solving the linear NPR model with the CPLEX algorithm. Table 5 gives the objective function values of solutions found by the CPLEX-based approach. In particular, results of the relationship score (last panel of the table) reveal that the solutions of the linear NPR model found by the CPLEX (NPR(L)-Cpx) and those by the tabu search (NPR(L)) are quite similar for instances with less complex constraints. With the maximum time limit set for each call to an optimizer, we could see that the solution returned by CPLEX is not necessarily optimal. Moreover, for relatively complex instances (pr11-12, [15][16][17][18], NPR(L)-Cpx failed to solve the problem and terminated the search without finding a solution within maximum time limit of 1800 seconds per subproblem. Table 6 and Fig 3 give the computation time of NPR(L)-Cpx versus NPR(L). It reveals that for instances with less complex constraints, the linear NPR model could also be solved using CPLEX. The computation time of CPLEX, however, increases significantly as the problem complexity increases. Consequently, for larger scale instances (in terms of constraints) the model can only be solved in a reasonable time with the proposed tabu search method.
Furthermore, the computation time of our proposed tabu search approach scales very well with the instance complexity, while CPLEX does not. In particular, calculating the

PLOS ONE
Nurse-patient relationship for multi-period home health care routing and scheduling problem improvement in the algorithm runtime (by taking the quotient of the runtimes NPR(L)-Cpx over NPR(L)), we obtained the ratio which ranges between 6x and 50x. This means our approach could achieve at least 6-fold runtime improvement over CPLEX. (For problems that CPLEX could not find a solution, this ratio is infinity). 5.3.3 Statistical tests on the relationship score obtained from the five models. Statistical tests have been carried out to see whether the obtained values of the relationship score found by the five model variants are significantly different. First, the obtained relationship scores of each model were aggregated in Table 7 by the mean-rank method. Note that when we calculated the mean rank, we only include the 14 problem instances for which the NPR(L)-Cpx  Table 7. Mean rank of the models based on the relationship scores, conducted on the 14 problem instances without N/A values. The higher the value, the better the model.

Model
Mean rank  Table 8 is significant, we further conducted a Wilcoxon signed-rank test to locate which pairs are significantly different at the 5% level of significance. The test results are summarized in Table 9. Indeed, the test found no significant difference between the linear and sigmoid NPR models. The test found also no significant difference between the linear models with CPLEX and the tabu search (N/A values removed from the calculation). NPR models are significantly better than the Basic and CC models. Finally, the Basic and CC models are not significantly different.

Observations of effects of the nurse-patient relationship score
Besides achieving a good relationship score, this section aims at illustrating the indirect effects of the NPR models through observations. While other objectives, namely, total distance and total preference score also influence the optimization results of the NPR models, as the models give priority to the relationship score, we present two easy-to-understand solutions to highlight the NPR effects on 1) minimizing the number of care workers and 2) promoting consecutive assignments. Figs 4 and 5 give a bar graph comparing the relationship score and an assignment table of a representative patient. Each figure consists of two parts: a bar graph and its corresponding table. The vertical axis represents the relationship score of the patient with the assigned care worker, and the horizontal axis represents the day index. The schedule table in the figure shows assignments where the column corresponds to day in the planning horizon, and the three rows give the care worker ID assigned to the patient. Note that if there is no request from a patient on any given day, then we do not display the value of the relationship score in the bar graph, and the assignment is marked by a dash (-). Recall that the CC model gives priority to minimize the number of care workers visiting a patient. While the priority of the NPR models is to maximize the nurse-patient relationship score, our results reveal that to achieve that goal the model simultaneously minimizes the number of care workers. The results showed that the relationship score of NPR on each day is highest. In addition, the assignment pattern indicates that NPR models have a tendency to assign care workers on consecutive days, while the other two models do not.
We see also that CC and NPR models both assign only Care workers 2 and 19 to the patient (despite being assigned in different order). Thus, the total number of different care workers for Patient 18 of CC and NPR models are both equal two care workers. Results on Day 19 however reveal different assignments made by the two models-the CC model assigns Care worker 2 because of a higher preference score (the preference score for Care worker 2 with Patient 18 is p 2 18 ¼ 0:941 while the preference score for Care worker 19 with Patient 18 is p 19 18 ¼ 0:633). This leads to a lower relationship score of CC when compared to NPR on Day 19 as shown in the graph. For the NPR models, they keep assigning Care worker 19 to the patient because Care worker 19 has been assigned to this patient previously for four consecutive times, thus having higher relationship score. This clearly demonstrates the essence of NPR as a function promoting consecutive assignments which should be one of the key indicators of quality of care-as using the same care worker can lead to the efficacy of monitoring patient's disease progression and treatment.
Let us remark that these direct observations of the effects of NPR score were exemplified by a fixed patient, fixed problem instance over the entire horizon. We are not attempting to generalize these observations to all patients by any means.

Conclusion and future work
The continuity of care model proposed in the previous studies of the multi-period HHC attempts to reduce the number of different care workers assigned to each patient without taking into account the consecutive assignments of the same care worker. Because of this, the CC model does not promote the consecutive assignments of the same care worker to meet with the patients. To systematically improve the CC model, this article presents a novel nursepatient relationship model. The dynamic nurse-patient relationship score captured by the NPR models demonstrate the trust a patient has for the care worker when they regularly meet. The model focuses on enhancing patient satisfaction by assigning a care worker with highest nurse-patient relationship score to the patient. To the best of our knowledge, there was no existing research work that examines a dynamic relationship function within the framework of HHC.
Extensive experimentation has been carried out. Twenty instances of the 28-day HHC problem are solved using the proposed heuristic algorithm combining a greedy algorithm and tabu search. Furthermore, to test the efficiency of the proposed approach, CPLEX has been used as an alternative solver to solve subproblems with less complex constraints. However, numerical results showed that CPLEX could not handle several problem instances with complex constraints, terminating the search without finding a solution. Tabu search, on the other hand, coped well with all problem instances. In addition, tabu search required almost the same computation time across all problem instances. The overall results indicated, in addition, that the NPR models improved care worker utilization, achieving highest relationship score for all instances. Moreover, not only the NPR models attempted to reduce the number of care workers assigned to a patient, we found also that the model promoted consecutive assignments arranging the same care worker to meet with the patient recently visited.
While our findings are important and promising, as a future work, there is still room for improving the sigmoid model to create a scoring system to reflect, capture closer to the realworld nurse-patient relationship, and investigate the consequences of applying the nursepatient relationship to the problem with other objective functions such as workload balance, number of tasks, fairness, etc. In addition, when dealing with uncertain scheduling problems, it will be interesting to consider the robust optimization for HHC planning as unplanned absenteeism among care workers might present especially during the pandemic due to the risk of COVID-19 infection.
In terms of limitations, although we have showed that the tabu search method could achieve much faster convergence than CPLEX when solving the NPR model, Table 6 reveals a drawback of the proposed NPR model compared to the CC model. It is evident from the table that computational complexity of the NPR models is a major bottleneck of practical implementation when trying to extend the model to a longer-horizon multi-period model. Another limitation is the fact that the trust level in the proposed NPR model only relies on the frequency of visits by a particular care worker. In reality, however, there are many scenarios that could contribute to relationship breakdown, e.g. work stress and burn out, nurse unfit to perform a job for a health condition, no call-no show at work, and several other possible causes. As these negative factors can cause a decrement in trust, modeling a more flexible relationship function that can incorporate other factors leading to the creation and/or breakdown of trust will be interesting avenues for future research.